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Local increases in the mean of a random field are detected (con- 
servatively) by thresholding a field of test statistics at a level u chosen 
to control the tail probability or p- value of its maximum. This p- value 
is approximated by the expected Euler characteristic (EC) of the ex- 
cursion set of the test statistic field above it, denoted E(^(j4„). Under 
isotropy, one can use the expansion E(/j(A„) = '^^Vkpk{u), where 
Vfe is an intrinsic volume of the parameter space and pk is an EC 
density of the field. EC densities are available for a number of pro- 
cesses, mainly those constructed from (multivariate) Gaussian fields 
via smooth functions. Using saddlepoint methods, we derive an ex- 
pansion for pk{u) for fields which are only approximately Gaussian, 
but for which higher-order cumulants are available. We focus on linear 
combinations of n independent non-Gaussian fields, whence a Central 
Limit theorem is in force. The threshold u is allowed to grow with 
the sample size n, in which case our expression has a smaller relative 
asymptotic error than the Gaussian EC density. Several illustrative 
examples including an application to "bubbles" data accompany the 
theory. 

1. Introduction. Data which can be modeled as realizations of a smooth 
random field are increasingly abundant, owing to advances in scanning tech- 
nology and computer storage. Traditionally, such data have arisen in the 
fields of oceanography, astronomy and medical imaging (see [40] for an easy 
introduction). More recently (e.g., [16] or the January 2005 front cover of 
Nature [3]), the "bubbles" paradigm in the cognitive sciences has provided 
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new and interesting applications. The difficult inference step is most often 
assessing the significance of the test statistic field arising from such data, as 
the multiple comparisons problem precludes "pixel-wise" significance test- 
ing. Commonly, the supremum of the statistic image is compared to an 
estimate of its upper pth quantile under the global null hypothesis of a zero 
mean function (equivalently, the statistic field is thresholded at this quantile 
to identify regions of activation) . A leading approach to approximating these 
quantiles uses the expected Euler characteristic (EC) of the excursion of the 
field above a threshold ([2, 38]), as described below. However, a common 
thread in these analyses is the assumption that the random field data are 
Gaussian; in some applications this assumption is known to be flawed. 

Practitioners of voxel-based morphometry in neuroimaging have high- 
lighted the problem [27, 37], but have stopped short of proposing a theo- 
retical solution. In the context of galaxy density random fields, Matsubara 
[21, 22] obtained approximations to the expected EC in three dimensions 
using Edgeworth series (though he never made the jump to a saddlepoint 
expansion). Catelan and coauthors [8] had previously used a similar ap- 
proach (also in an astrophysics context) to get at the mean cluster size of 
excursions for non-Gaussian fields. Rabinowitz and Siegmund [25] used sad- 
dlepoint methods to approximate the tail probability for the supremum of a 
particular non-Gaussian field, a smoothed Poisson point process. Although 
they supplied a first-order version of the expected EC and published heuris- 
tic arguments, their work originated some of the central lines of reasoning of 
the current article (indeed, since the paper is cited often we henceforth refer 
to it as RS). We shall formulate a general approximation to the expected 
EC for asymptotically Gaussian fields with known (decaying) cumulants, 
and compare its relative error to that under a full Gaussian assumption. 
In particular, our approximation is found to outperform the Gaussian ex- 
pected EC in a particular "large threshold" asymptotic setting. Throughout, 
we underscore the connection between our results and those of RS. 

For concreteness, we consider processes of the form 



where the Wj's are i.i.d. smooth random fields on T C with zero mean, 
constant variance <t^ and finite higher cumulants. Such a Z{t) is a special 
case of what we call a Central Limit random field (CLRF), and usually 
should be thought of as a test statistic derived from data. We call Au = 
Au{Z,T) = {t £T; Z{t) > u} the excursion set of Z above u. We denote the 
EC, an integer- valued set functional, by ^{■). For a working definition of the 
EC, see [34], and for a more rigorous account the reader is referred to Part 
II of Adler and Taylor [2]. As the latter volume is the second major source 
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Fig. 1. Galaxies are not uniformly distributed in the universe, but tend to congregate in 
clusters, strings or even sheets. Cosmologists have used the EC of regions of high galaxy 
density to characterize this large-scale structure ([18] used data from the first release of 
the Sloan Digital Sky Survey SDSS). In 2D the EC counts # "blobs" — # "holes" in a set. 
We examine luminous red galaxies (LRC's — the most distant objects in the SDSS) from 
the latest (5th) data release [24]- In (a) we show LRC's from a 22.5° square patch at 
850 ± 50 megaparsecs (Mpc) with roughly the same galaxy density as in [18] (our sample 
is part of a sphere, rather than cone-shaped, so we avoid the luminosity correction). Panel 
(b) displays the same number of "galaxies" simulated uniformly according to a Poisson 
process. In (c)-(d) the LRC's are smoothed with a Caussian kernel of s.d. 10 Mpc, and 
standardized to have zero mean and unit variance under the Poisson model. 



of methodology used in this article, it will be referenced as AT. Very loosely 
(ignoring what happens on the boundary), in 2D the EC counts 7^ "blobs" — 
# "holes" in a set, while in 3D it counts # "blobs" - # "tunnels" + # "hollows" 
[40]. 

In some applications the expected EC of A function of u is itself 

of central interest. For example, cosmologists compare the observed EC (or 
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Fig. 1. (Continued.) In (e) we compare the observed EC of each image at a continuum 
of thresholds with the expected EC under the Poisson model computed using a Gaussian 
approximation and the "tilted" correction derived in Sections 4-5. The observed EC for the 
SDSS image deviates from the others, reinforcing the nonuniformity of LRG's. The Gaus- 
sian and tilted expected EC's differ primarily in the tails; in the range where the expected 
EC approximates a p-value (bottom-right zoom) this correction is dramatic. In particular, 
the tilted curve majorizes the Gaussian one since the LRG distribution has positive cu- 
mulants. The observed EC for the (asymptotically Gaussian) Poisson data conforms well 
with both curves. Note that all curves approach 1 as u ^ —oo and as u ^ oo (top-left 
zoom). 

1 — EC, sometimes called the genus) from galaxy sm'vey data with the ex- 
pected EC from various hypothesized models of the large-scale structure of 
the universe (see [18] and the illustration in Figure 1). In other applica- 
tions, including neuroimaging and the bubbles problem treated in Section 
8, interest centers around using this quantity as an approximation to the 
p-value: 

(1.2) Ev9(y4„) «p|supZ(i) >iil 

for large u. The assumption that (1.2) — the so-called "EC heuristic" — is 
adequate for statistical inference is no longer a mere heuristic for Gaussian 
fields; see [33] or AT, Chapter 14. See also Section 4.4 of the present paper. 

To motivate the study of what are known as EC densities, we have the 
following well-established result. 

Theorem 1 ([39], AT). Let Z{t) he any isotropic random field on T d 
M^, a locally convex finite union of convex bodies. Then the expected EC of 
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its excursion set at the level n G M can be written 

N 

(1.3) Eip{Au) = Y.Vk{T)pk{u), 

k=0 

where Vfc(T) is the kth intrinsic volume of the search region and pk{u) is a 
function called the kth EC density of the field. 

Proof. The expected EC E(p{Au{Z,-)) is a real-valued additive func- 
tional acting on T; Eip{Au{Z,T U V)) = Ecp{Au{Z,T) U Au{Z,V)) = 
Kip{Au{Z, T)) +Eip{Au{Z, V)) - Ecp{Au{Z, T) n Au{Z, V)) = Eip{Au{Z, T)) + 
Kip{Au{Z,V)) — Eip{Au{Z,T ri V)). If the process is isotropic, then 
Eip{Au{Z,-)) is also translation- and rotation-invariant. Finally, the local 
convexity of T guarantees that the functional is continuous (e.g., in the 
Hausdorff metric). The result then follows from Hadwiger's seminal theo- 
rem [17, 28], which states that any such functional is in the linear span of 
the intrinsic volumes, and hence admits an expansion like (1.3). □ 

Remark. If the field is nonisotropic but derived from a (multivariate) 
Gaussian field by a smooth function, then on sufficiently regular parameter 
spaces (1.3) holds with Vk{T) replaced by the kth Lipschitz-Killing curva- 
ture Ck{T). The latter differential geometric functionals depend both on the 
space T and on the covariance structure of the field itself; see AT Chapter 
10. 

1.1. Intrinsic volumes. The intrinsic volumes of a set T are a general- 
ization of its volume to lower dimensional measures. There are several ways 
of describing these functionals; here we give an implicit definition, which is 
valid for any locally convex set T and sufficiently small radius r. Let | • | 
denote the Lebesgue measure, B{0,r) be the ball of radius r centred at 0, 
and T © -6(0, r) = {x + y.x £T,y G B{0, r)} be the tube of radius r around 
T. Write oJk = \Bk{0, 1)| = n''/'^ /T{k/2 + 1) for the Lebes gue measure of the 
unit bah in M^. Then 

N 

(1.4) \TeBiO,r)\ = J2^N-kr''~^VkiT). 

k=0 

Thus, the intrinsic volumes are related to the coefficients in the Steiner- 
Weyl volume of tubes expansion. For example, in M.^ we have Vo(T') = ^{T)^ 
Vi(r) = 2 X caliper diameter (T), V2(T) = 1/2 x surface area(T), and V3(r) = 
volume (T). Remarkably, EC densities for all "Gaussian-related" processes 
were shown to have a similar characterization in [31]. 
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1.2. Computing EC densities. For a matrix A or vector a we introduce 
the notation j4|fc = (j4jj)i<jj<fc and a\k = (ai)i<j<fe- When k = this de- 
notes, by convention, an "empty" matrix or vector — that is, a term or factor 
which disappears from the expression. We use Ip (or just / when the dimen- 
sion is unambiguous) to denote the pxp identity matrix. When apphed to a 
random field, "'" represents spatial differentiation with respect to t, and we 
shall refer to the individual components of Z and Z, respectively, by Zi and 
Zij . There are several useful formulations for computing the EC densities of 
stationary fields [39] . In the derivation of Theorem 4 we use 

(1.5) pk{u) =E{Zk^[z^:>o]<^et{-Z\k^l)\Z = u,Z\k^i = 0}/o,fc-i('u, 0) 

for k = I, . . . , N , where /o,fe-i(-) denotes the joint density of {Z, Z\^_^. In- 
deed, this was the original method used by Adler to compute the expected 
EC, and corresponds to using the A^th coordinate function (ti, • • • , ^Af) *■ 
t^ (instead of the field itself) as the function in Morse's Theorem (see [1], 
Chapter 4). The zeroth EC density is just the univariate tail probability, 

Pq{u) = P{Z(t) > u}. It has a clear interpretation as the lowest order ap- 
proximation to the p- value (1-2). 

1.3. A glimpse ahead. EC densities are available in diverse parameter 
spaces for a number of processes [6, 7, 30, 32, 35, 38], mainly those de- 
rived from (multivariate) Gaussian fields via smooth functions. Let Hm{x) = 

ml X]j=o^^ 2^j\{m-2jy. mth Hermite polynomial. For a mean zero Gaus- 

sian field with variance o"^, the densities take the form (cf. AT, Chapter 11) 

(1.6) p^(.)^(2.)-(^+i)/V-^exp{-^}i7,_,(^). 

(Note that 7 is used to denote the Gaussian measure, following the con- 
vention of AT.) Our main result (Theorem 4) will be to show that for a 
stationary CLRF like (1.1), and assuming that the threshold u is of the or- 
der for some l/4<a<l/2, modifying only the exponential portion 
of (1.6) leads to an approximation to pk{u) whose relative error behaves like 
^i/2-2«^ Adding factor related to the mixed skewness terms E[Z^Zj] can 
yield an additional improvement. Moreover, we shall show that this error is 
asymptotically strictly better than that of (1.6), exponentially so over the 
range 1/4 < a < 1/3. 

The remainder of the paper is organized as follows. In Section 2 we use 
saddlepoint methods to approximate the joint density of {Z,Z' ,Z') (suit- 
ably transformed), leading to a mixed tilted-direct Edgeworth expansion. In 
Section 3 we consider the growth rate of u and compute certain moments 
of the tilted distribution which shall be used in subsequent derivations. We 
derive the approximation to Pkiu), by methods reminiscent of [1] and [38], 
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in Section 4. Section 5 offers an analysis of the relative error of the tilted 
densities compared to that of a standard Gaussian approximation, and we 
comment on the contribution of RS in Section 6. In Section 7 we illustrate 
the method by treating a field — the Xn process — whose densities are known 
exactly, and Section 8 applies the theory to the test statistic from a bub- 
bles experiment. We also present a simulation under the latter framework 
which examines the accuracy of our corrected version of (1.3) as a p-value 
approximation. In the interest of brevity, we have omitted several proofs and 
technical developments, all of which can be found in an online version of the 
present article (URL), or in the doctoral thesis [9]. 

2. A mixed expansion. In order to avoid having EC densities which de- 
pend on t in Theorem 4 below, we shall assume that Z is strictly station- 
ary. For practical purposes [i.e., to apply Theorem 1 in the computation 
of E</j(74„)], the stronger requirement that Z be strictly isotropic shall be 
in force. In that case, purely for convenience, Z will be assumed to have 
unit second spectral moment, namely M[ZiZj] =5ij. The latter assumption 
can easily be relaxed to allow Var[Z(t)] = XI, which merely modifies (1.3) by 
multiplying the kth term by A'^/^. Strict isotropy can, on very simple param- 
eter spaces, be weakened to stationarity [with a similar correction to (1.3)]; 
see Section 4.3 or AT, Chapter 11. If > 0, the notation x„ = 0{n~^) will 
indicate that for some < C < cxd and all sufficiently large n, \xn\ < Cn~^. 
A similar notation is used for matrices, with the understanding that the re- 
lation holds componentwise. The notation a„ ~ bn signifies that an/bn 1 
as n — > oo. 

Fix k>2 and put 

(2.1) M = Mfc = aZ|fc„i + -4_i. 

a 

The most important property of this shifted (A: — l)x(A; — 1) Hessian is 

that it is uncorrelated with the field: K{MijZ} = 0. Let Xijim = ^{ZijZim} 
[a symmetric, 0(1) function] denote a fourth-order spectral moment of the 
process. It is easily shown that 

We shall write M to denote both the matrix and the corresponding random 
vector vec(M), with the context eliminating confusion (similarly for Z). In 

view of this, let Z = = {Z, Z\'f,, M')' have density /(z), to which we shall 
apply the asymptotic expansion. Define the exponential family 



(2.2) fe{z) = fe{z, z, m) = exp {Oz - K{9)]f{z, i, m). 



N. CHAMANDY, K. J. WORSLEY, J. TAYLOR AND F. GOSSELIN 



where 6 = {Oo, . . .,Od^^i) G = k{k + l)/2 + 1. The function K{-) is of 
course the cumulant generating function (cgf) of Z, 

oo 1 

(2.3) i^w = E E -kA[mo\ 

\"/\=U 

71 7d 

where if 7 G we define [Y]^ = ^^T^, • ■ • ,^7^^, (^i, • • • , e^V = 
I7I =J2i=i7i^ and 7! = Ui=i7i-- We have written k^, for a z^th 
order cumulant. 

Tilting in the first argument and inverting (2.2) yields 

(2.4) f{z,z,m) =exp{K(6'o,0') - 9oz}f(^g^^^o'){z,z,m). 

The relationship (2.4) holds for all values of Oq, in particular, that value 
^0 = which makes the Edgeworth expansion to /(g^ q/)(z) most accurate 
at the point z = (ti,0',0')' — namely, that for which 

5^(^0,0') 



(2.5) ^ = E,..oM^ 



holds. It is trivial to see that this requirement is satisfied if and only if 6^ 
is the formal maximum likelihood estimate of 9q under the model /(g^, o')(^) 
when Z = It is observed. 

Finally, we let fi/ and denote the mean and covariance matrix under 
the density f^g o')^'^^' "^^^^ (cf- [4], page 187) we have 

^V=(v,....,K(^)|,^(,^,,J ^/=^^'m\e^0^,o'y 

where V2-d^ has been used to denote {8/ 862-, ■ ■ ■ , 8/89df,)' ■ The approxima- 
tion to be exploited will eventually take the form 

f{z,z,m) =exp{K{0u,O') - (9„z}0d^((z,i',m')';Ai/,S/) 
(2.6) X {1 + Qk,3iiz, i',m'y - fi/)} + ek,n, 

= fk{z, z, m){l + Qk,z{{z, i', m')' - /x/)} + efc,„, 

where Qk,2, is a degree 3 polynomial containing mixed cumulants of (Z, Z\k,M) 
of order 2 and 3, while ek^n denotes terms of smaller asymptotic order. It 
turns out that Qk,^ has coefficients (including constant term) which range 
up to 0(n~^/^). Expansion (2.6) is called "mixed tilted-direct" because ex- 
ponential tilting is applied to the first argument but not to the others. 
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3. Preliminary calculations. We argued heuristically that the expected 
EC will provide a good approximation to the p-value (1.2) for large thresh- 
olds u] we therefore allow the possibility that u = Un ^ oo as n^oo. Let 
u = 0(n^/^~") for some < a < 1/2. Prom (2.5) we have 



(3.1) 



where Kz{-) denotes the marginal cgf of Z. When is a univariate ar- 
gument, we shall use ""' or a superscript [e.g., "(^)"] to denote (multiple) 
differentiation with respect to 0. Under some additional assumptions, (3.1) 
is a proper asymptotic expansion, in the sense that the error in the remain- 
der is of smaller order than the previous term. We shall use the following 
lemma repeatedly, and often implicitly, in proving our main results. 

Lemma 2. Let Z = n^^^'^J^i^i ^^^^ centred, i.i.d. random vari- 
ables with variance such that E,\Wi\'^~^^ < oo for some s > 1. Suppose 

that K^'^^^ is continuous on a closed interval [—eg, eg] with > 0. Let 
\9\ < cn^^'^~^ for some finite c and /? > 0. Then there exists a constant c{s) 
depending only on s such that, for all large n, 

q2 qs~1 



(3.2) 



K'z{e)-a'e-k^{Z)- kg{Z)- 



< c(s)ni/2-«/3. 



2 {s-iy. 

In particular, if u ^ cn^/^""^ for some < a < 1/2 and c > 0, with Z as 
described above and s = 2, then 

(3.3) eu = eu,n = ^[l + 0{n~'')\ 

as n—> CO. If a = 1/2, then u is (asymptotically) constant and 6^ = u/a"^ + 
0(n-V2). 

3.1. The tilted moments. Let us now be more explicit about the "tilted" 
parameters /// and E/, computed by differentiating (2.3) with respect to the 

appropriate argument(s) and then setting 6 = {9u,0')- To simplify notation, 
we partition as follows: 

We alert the reader that the dot notation in the above parameters does not 
indicate differentiation, but was designed to emphasize that portion of the 
process Y{t) with which a particular mean or covariance is associated. In S/, 
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the diagonal entries in the partition are just the tilted variances, respectively, 
of Z, Z\k and M. The off-diagonal entries are the corresponding covariances, 
and so CT is /c X 1, (T is k{k — l)/2 x 1, and S is k{k — l)/2 x k. We have, for 
example, 

OO -1 

i>=3 ^ '' 

and many more such expressions. When a > 1/4 it is easily deduced that 

(3.5) rl^T^ = a^ + 0{n-^), a, = 0(?i""), <t,,- = 0(n^-); 
S,,,7 = 0(n^"). 

3.2. Some tilted conditional moments. We shall also need the following 
conditional moments (assuming a > 1/4) under the tilted Normal density 

4^k{-) = 0dfe(s A*/) 5]/), the proofs of which are straightforward: 

(3.6) u = E^jZfclZ = u, Z\k-i = 0} = 0(ni/2-2°); 

(3.7) = Var^jZfclZ = n, = 0} = 1 + ©(n^"); 



(3.8) 



(3.9) 

where Cij = 0{n~°') and 



IE</.fc{Mjj|Z = n,Z|fc_i = 0,Zk = x} 
= /ii, + Q,x + 0(nV2-3a). 

E^jMijMi„|Z = u, = 0, Zfc = x} 



(3.10) 4,^ = 0(nV2-3a), 

4,^ = 0(n-2"). 

The punchline to (3.8) and (3.9) is that under the Gaussian approximation to 
the conjugate density /^^ and conditional on the event {Z = u, Z\].^i = 
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0, Zk = x}, we have M = A + r-x, where A is a centred Gaussian matrix with 
covariances satisfying (4.1) below, and H^,. is a smah, nonrandom (given x) 
perturbation. 

4. Euler characteristic densities. We shall now proceed to use the char- 
acterization (1.5) in order to derive our corrected EC density. 

4.1. The case a> ^. The following analogue of AT's Lemma 11.5.2 can 
be proved in a similar fashion if one is mindful of the various error terms. 
In what follows, | • | is the determinant. We also add a subscript n to all EC 
density expressions, since we are interested in their asymptotic behavior. 

Lemma 3. Let \ < a < ^, let < x be a 0(1) nonrandom scalar and 
£ a 0(1) symmetric function from {l,...,fc}^ into M. Suppose that A is a 
centered, (jointly) Gaussian k x k matrix whose covariances satisfy 

(4.1) E{AijAi^) =e{i,j,l,m) - 6ij6l„^ + Siji„i{x), 

where Sijimix) is a quadratic polynomial with coefficients as in (3.10). Let 
h = 0{n^^'^~'^) he a nonrandom scalar, and let = 0(n^/^~^°), r} = 0(n~"), 
and = H'' + xr} he nonrandom matrices. Then 

(4.2) E| A -hl + = {-lfHk{h){l + rfc(x)}, 

where r^ is a degree k polynomial whose largest (constant) coefficient is 
aQ = 0{n~") as oo. 

Theorem 4. Let Z he a strictly stationary CLRF as in (1.1), with co- 
variance function belonging to C"^(]R^) and E[ZZ'\ = I . Let u ~ cn^/^~" as 
n oo for some constants < c < co and \ < a < ^, and suppose that 6^ is 
a solution to 

(4.3) K'zi9) = u, 

where Kz{-) is the cgf of Z{t). Let = K'^iOu) he the tilted variance of Z{t) 
at 9u, and Iu{9) =u9 — Kz{9). Then subject to regularity conditions, the EC 
densities of Z for > 1 are given by 

(4.4) /Ofc,n('w) = Pk,n{u) X {1 + e„}, 
where 

(4.5) pk,n{u) = {2TT)-^^+^y\-^eM-Iu{0u)}Hk-i{ru9u) 
and as n—> oo the relative error satisfies e„ = 0(n^/^~^"). 
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Moreover, an approximation with relative error 0{n^^^^ "'^ ^"^) is given 

by 

(4.6) p,,„(n)x (l + E[z2z,]y||), 
and (if k>2) one with relative error (9(72™ax{-a,3/2-6a}-j 

(4.7) p,,n{u) X (l + K[z'z,]^l^^-^)exp{-^4T.nz'Z,f]. 



1=1 



Remark. The "regularity conditions" mentioned above are just those 
on the random field (Z, Z|'^_^) required to apply the Expectation Meta- 
Theorem (Theorem 11.1.1 in AT), those on the density f(z,z,m) assuring 
the correctness of the Edgeworth expansion, plus a mild assumption about 
the moments of {W, W, W). These conditions are academic in that they can 
rarely be checked in practice, but are considerably weaker than the Gaussian 
assumption, which is currently the only alternative. See the Appendix for 
further details. 

Proof of Theorem 4 (Leading term). For extra details the reader 
is referred to the online version of this article, where we rigorously bound 
the error terms and consider the contribution to the integral (1.5) of the 
polynomial Qk,3{')- It is seen there that the additional terms are of relative 
size at most 0(n~"), and are therefore negligible. The crux of the derivation, 
however, is in working out the leading term, which we do now. Plugging the 
saddlepoint approximation into (1.5) yields 

(4.8) pf, {u)= / xdeti )fk{u,0',x,m)dmdx 

= exp{-Iui9u)}{ 



a 



X 

(4.9) 



/■oo 

/ x(t)k{x\u,Qi') 

Jx=0 

X / deiim-u/aI)Mm\u,0',x)dn.dx 
X 4'o,k-i{u,0'), 

where Dj = j{j + l)/2 and the various (j)'s denote the obvious Gaussian 
densities derived from p/,T,/) by conditioning or marginalizing. Eval- 

uating the inner expectation via (3.8)-(3.9) and Lemma 3, after noting that 
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(4.10) 



Tu^u — u/a = 0{n °) has smaller order than rP , gives 
exp{-/„(a,)}a-('^-i)//fc_i(r,^„) 

xE{y(l + rfc_i(y))l[y>o]}<Ao,fc-i(^,0'), 

where Y r-. M{v,rf) as in (3.6) and (3.7). Recall that cr = r„(l + ©(n"")), 
and consider the last factor in (4.10). We have 

(4.11) </>o,fc-i(n,0') = (27r)-^-/2det(Ffc)"i/2exp{-iAlLiS|^^i/iU-i}, 



where 
(4.12) 



Vk 



So det(Vfc) = r^(l + 0{n ")), and the exponential portion of (4.11) behaves 
like 

exp{0(ni/2-2ay[j ^ 0{n-'^)]0{n^'^-^'')] = exp{0(ni-^")} 

= l + 0(n^"^°). 
It is also standard fare (e.g. [1], Lemma 5.3.3) that 



E{yi[y>o]} = ^exp 



27r 



(•- 



2'rf') ^ \^rj 
l-|i + 0(n--) 



'2-n 

+ 0(n3/2-6'^) 

1 _|_ ^ _|_ Q|^^max(-Q,3/2-6a)^ 



2 ^2^?? 



(4.14) =(27r)-i/2^0(ni/2-2a)_ 

Let pfc(2/) = yrk-i{y) and = (^ - i')/??- Write pk{y) = ELi aj-i?/', with 
the aj's at most 0(n~"). Then 



(4.15) 



3=0 



■i=max(l,j) 



j=0 
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where the largest coefficient cj in (4.15) can be loosely bounded by 0{n'^ 
The absolute moments of Yq are all finite, so that 

k 

\Ep,iY)l[Y>o]\ < E |cj|IE|loP = 0(n-"), 

j=0 

whereby 

E{y(l + rfc_i(y))l[y>o]} = (2vr)-^/2 + ©(n^/^-^-). 
Putting the pieces together, (4.10) simplifies to 

e"^"(^")r„-^(2vr)-('=+^)/2F,_i(r„^„) 

X (1 + 0(ni/2-2"))(i + 0(n-"))(l + 0(ni-^")) 

= PMWx(l+0(ni/'-'")). 
Toward (4.6), consider one additional term in (4.13)-(4.14): 

(27r)-i/2 ^ + 0(n"^^^^-"'3/2-W) 

= (27r)-i/2 + ^^/2 + 0(n'^^^^-"'3/2-6a}) 

= (27r)-l/2 ^ ^^(^^ Zfc)^2/4 + 0(„max{-a,3/2-6a}) 

32 



(27r)-i/2|i + E[Z2Zfc]y^^|(l + 0(n'"^^^-"'3/2-6"})). 



As for (4.7), we have 



-^Alfc-iSlfc^i/ilfe-i = + 0{n "))/iU-i 

= -i E (E[z2z,]^2/2 + 0(ni/2-3"))2 

i<k-l 



+ AlLiO(n-")Ak.-i 



Finally, note that 

max{l/2 — 3a, 1 — 5a} < —a < max{— a, 3/2 — 6a} 

< max{— a, 1 — 4a}. □ 

It is worth emphasizing that if the mixed (third-order) cumulants ¥,[Z'^Zi] 
of Z with its first derivatives vanish, then the approximation is improved 
regardless of the behavior of the marginal cumulants of Z. Note that if 1/3 < 
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a < 1/2, then (4.7) is no better than (4.6), since —a > 1 — 4a. An important 
special case of Theorem 4 is when a = 1/2, so that u is a constant. In that 
"small-threshold" setting each of the approximations (4.5)-(4.7) admits a 
relative error of the order 

4.2. What about a < j? Consider the case a < j, in which the threshold 
u grows somewhat more quickly relative to n. Here we expect the EC heuris- 
tic to perform well, and the Gaussian approximation to perform poorly. But 
what of the tilting procedure? From (3.5) we have that when a < 1/4, the 
tilted mean derivatives also grow with n: fii, jlij ^ oo. Lemma 3 no longer 
holds in this scenario since the "remainder" terms now dominate (4.1) and, 
thus, the special covariance structure which leads to Hermite polynomials is 
lost (the problem is exacerbated if a < 0). 

A reasonable conjecture in this setting is that the tilted EC densities will 
provide a much better approximation to the true ones than will those derived 
from a Gaussian assumption (see the next section). Less obvious is whether 
the tilted densities themselves have an acceptable relative error — indeed. 
Theorem 4 seems to suggest the contrary. While we have yet to work out 
the details in general, empirical evidence in at least two cases (see Sections 
6 and 7) supports using the tilted densities. 

4.3. Brief comment on Rabinowitz and Siegmund. In the language and 
notation used here, the Rabinowitz-Siegmund heuristic can be phrased 

(4.16) p|supZ(t) > u| w 1 -exp{-E{A4(Z)}} wE{M„(Z)}, 

where Mu{Z) is the number of local maxima of the field Z above u. In RS, 
the authors approximated this expectation by 



(4.17) E{M,(Z)} = exp{-/„(0J} 



(2vr)(fc+i)/2r„ « 



where A = Var^ [Z] is the tilted variance of the full derivative process (matrix 
of second spectral moments), having removed the local isotropy condition. 
Note that 0^~^/tu is the leading term in the polynomial Hn-i{tuOu), 
so that (modulo the spectral determinant) (4.17) is equivalent to the leading 
term of our approximation. Naively extending (4.17) for, say, a stationary 
process on a rectangle, one might obtain 

N 

(4.18) Eip{Au) = J2 E \J\k\h\^^%,niu)^ 

where now Aj denotes the tilted variance of the derivatives restricted to 
facet J of the rectangle. Ok is the set of A; -dimensional facets touching the 
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origin (cf. AT, Chapter 11), and | • \k is a /c-dimensional Lebesgue measure. 
(See Section 5.2 for a definition of po,n-) For a general isotropic process, 
where A = XI, (4.18) would become 

N 

(4.19) E(^(A„) = ^ A^/2Vfc(r)pfc,„(n). 

It is readily seen that the relative error between the approximations (4.18)- 

(4.19) and our own versions (with the untilted A or Aj replacing A or Aj) 
is of order 0(n-"). Indeed, A^-; = E{Zj{t)Zi{t)} = 0(1), and 

oo 

(4.20) Xji = ^ K{Z, ...,Z, Z„Zi)j-^^ = A,,(l + 0(n-")). 

Hence, |A|^/^ = |A|^/^(1 + 0(n~°)), and similarly for Aj. This implies a 
0(n~") discrepancy for each A; > 1 term in (4.18) and (4.19). Since —a < 
1/2 — 2a, it is unclear from the theory whether one should use the tilted 
or untilted second spectral moment(s) in computing Eip{Au)- Similarly, it 
is seen in Section 5 and the proof of Theorem 4 that using the Gaussian 
argument u/a in the place of TuOu in the Hermite polynomials leaves the 
relative error unchanged from (4.4). Intuition suggests that TuOu may give 
more accurate results in finite samples. A geometric argument also supports 
its use; see [11]. The A; = term is unaffected by these changes, and can be 
improved by the methods to be discussed in Section 5.2. 

4.4. Some remarks about validity. It has been shown that for a Gaussian 
process the relative error of the expected EC as an approximation to the 
p- value is exponential as w — > oo, in the sense that the first term of their 
difference is exponentially smaller than the zeroth term in expansion (1.3), 
which behaves like u-^e-^'^/^^aa) [33]^ r^^^ authors of [34] (and an anony- 
mous reviewer) rightly pointed out that this result is, for the time being, a 
purely Gaussian one: the exponential error likely does not hold even for a t 
random field. Nevertheless, much of the development in [33] makes no use 
of Gaussianity. We argue here that for CLRF's, the expected EC is expo- 
nentially close to the expected number of local maxima above u, which is 
obviously an upper bound for the p-value. 

Assuming that the parameter space is a manifold without boundary, one 
can write the difference (before taking expectations) as 

(4.21) \ip{Au) -M^\<#{teT: Z{t) > u, Z{t) = 0, Z{t) ^ 0} 

(4.22) = G r : Z{t) > u, Z{t) = 0, {Z{t) + Z{t)I) ^ Z{t)I} 

(4.23) < #{t G T : Z{t) > u, Z{t) = 0, Ci{Z{t) + Z{t)I) > u}. 
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where Ci(') extracts the largest eigenvalue. We have again assumed VarZ(t) = 
/, purely for convenience. It is not hard to argue (e.g., using AT's Expecta- 
tion Meta-Theorem in combination with an Edgeworth expansion) that for 
a stationary CLRF the expectation of (4.23) behaves to first order like 

|2-|g-«V{2<x2)g-«V{2<x?)^ 

where a1 = supy^n^]^ Var[?;'(Z + ZI)v\ is the so-called critical variance. The 

key fact is that Z + ZI^ Z and Z are uncorrelated, and indeed almost inde- 
pendent. When T has a boundary, (4.21) must account for points where 
intersects dT] the argument is more finicky but not fundamentally different. 
In the interest of brevity, we leave the details to a future article. 

Since p < EM^ ~ Eip{Au), clearly using the expected EC is conservative 
(at least asymptotically). Determining conditions under which the p- value 
approximation is sharp will require modification of the purely Gaussian por- 
tions of [33]. Numerical evidence from a simulation (see Section 7) is reas- 
suring. 

5. Comparison with the Gaussian approximation. We shall use the no- 
tation 

(5.1) en[cn,dn] = (3 ^ c„ = d„ X ( 1 + O(n^) ) 

as n — > oo to describe relative error for (3 <0. We say that the error is sharp, 
and write e„[cn,(in] = /3, if (5.1) holds, but fails if /3 is replaced by /? — e 
for any e > 0. It is important to determine whether there are values of a 
for which our derived expression provides a better approximation to Pk,n{u) 
than does the Gaussian EC density. We examine only the approximation 
Pk,n{u)', similar considerations for the refined formulae (4.6) and (4.7) are 
straightforward and left to the reader. 

5.1. EC densities with A; > 1. In Theorem 4 we derived a "tilted" order 
0(re^/^~^") (in relative error) approximation to Pk,n{u) for A; > 1, denoted by 
Pk,n{u) and given in (4.5). Comparing this approximation to the Gaussian 
EC density p\ = p\n from (1.6), we have (under the regularity conditions of 
Theorem 4) the following: 

Corollary 5. Let u ^ cn some 0<c<oo. If ^ < a < and 

A;3(Z) 7^ 0, then for all k>l, 

(5.2) en[pl^n{u),pk,n{u)] = l-3a. 
In particular, 

(5.3) en[/Ofc,„(ii),Pfc,n(")] - e„[pfc,„(ii),pfc,„(u)] > ^ - a > 0. 



18 N. CHAMANDY, K. J. WORSLEY, J. TAYLOR AND F. GOSSELIN 



Moreover, if ^ < a < and k3{Z) / 0, then the error of the Gaussian den- 
sities is exponential, in the sense that 

(5.4) Pk,ni^) = Pk,n{u) X exp(6„) 

for A; > 1 and some sequence with |6„| ^ oo. 

Remark. In the proof of Corollary 5, it is easily seen that the sequence 
referred to in (5.4) can be written 6„ = k'i{Z)9^/Q + 0{n^~^'^), where the 
skewness may be positive or negative. 

5.2. The k = case. Until now we have only considered EC densities for 
k>l. The zeroth density is generally defined by 

(5.5) po,niu)=F{Z{t)>u}, 

the univariate tail probability above the level u. Since nothing in the deriva- 
tions of Theorem 4 or Corollary 5 explicitly requires u > 0, the results men- 
tioned thus far are valid if we weaken the threshold condition to u ~ cn^^^"", 
|c| < oo. [Their use is limited, however, by the facts that (a) when \u\ is small 
the Gaussian approximation is adequate and (b) when u is large and negative 
the zeroth density dwarfs the others.] With the exception of the Lugananni- 
Rice formula (5.9), the discussion in this subsection does assume n > 0, since 
it pertains to saddlepoint-type approximations to (5.5). If one wishes to ap- 
proximate po,n(^) for n < (as is the case in Figure 1), then the expressions 
can be modified by making use of the symmetry in their derivations and of 
the fact that po,n(^^) = 1 — ^{Z (t) < u} . Such generalizations are left to the 
reader. It should be noted that when a p- value approximation is sought and 
hence u is large, the error in po,niu) typically has a negligible impact com- 
pared to those in the higher EC densities. Approximating tail probabilities 
via saddlepoint methods is far from a new problem; the principal purpose 
of this subsection is merely to enumerate some useful formulae. 

Let ^(•) = 1 — $(•) denote the standard Normal survival function. Given 
Theorem 4, one who is inspired by AT's definition H-i{x) = \/2TT^{x)e^ 1"^ 
might naively estimate (5.5) by 

po,n(n) = exp{-/„(0"„)}(27r)-V2^_^(^^^"j 
= exp{r2^1/2 - /„(^„)}^(tJJ. 

Perhaps surprisingly, this is what Robinson [26] called the "first" saddlepoint 
approximation, and turns out to be a reasonable estimate if u > 0. The zeroth 
Gaussian density is 

PlM ^ (1 + O(n^--i)). 
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Under typical Central Limit Theorem conditions (e.g., [15], Theorem 2.3), 
one has to first order 

,2 



(5.6) 



u 



lU{u/a), 



/50,n(n) = PoV«(^) + ( "^2 

which implies (sharply, assuming nonzero skewness) that 

Daniels [14] gave two differently-derived formulae. The first can be written 
Po,n(^i) = exp{r^^^/2 - /„((9„)} 

(5.7) 



^(.A)(l-™« 



K 



(3) I 



(rX - 1) 



6 



and corresponds to the so-called "second" saddlepoint approximation [26]. 
Let /z(-) denote the marginal density of Z{t). Upon using the large x esti- 
mate "^(x) = (f){x)/x X (1 -|-0(x~^)), and under the regularity condition (cf. 
[26]), 



(5.8) sup 

z 

(5.7) leads to 



z — u 



PO,n{u) =P0,n{u)il+O{n °)). 

The second approach is the Lugannani-Rice formula [20] 

1 1 



PO,n{'^) = *(^«) +<P{^u) 



(5.9) 



Cu 



+ 



1 (Kf{e^) 5K§\e^f 



K 



(3), 



2C^ 
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+ 1 I + 



where = sgn[0„] y 2/,j(0m) and Qu = Tu^u- Note that this approximation is 
valid for large negative as well as positive n, and consequently, the first two 
terms in (5.9) were used in the astrophysics example of Figure 1. One can 
reduce (5.9) to 

Po,n{u) = ^{i0u){l + O{n-'')) 



(5.10) 



*(sgn[eJV2^«(^«))(l + 0(n"")) 



^^o,n(^)(l + 0(n-"))- 
Filling in a few details gives the following: 
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Proposition 6. Let \<a<^ and suppose that (5.8) holds. Then 

(-n[PQ,n{u),PQ,n{u)] = £„ [po,n H , PO,n («)] = -O. 

Also, if the field Z is as in (1.1), ks{Z) ^ and \ <ol<\, then 

en[/5o,„N,PO,n(w)] = 1 - 3a. 

In particular, 

en.[Po,nN'PO,nW] > e„, [po,n N , PO,n Va G [1/3,1/2). 

When a < 1/3 the error in the zeroth Gaussian density diverges to infinity. 

So once again the tilted densities perform better than the Gaussian one. 
However, po^n and /3o,n can both be improved upon if one is prepared to 
numericahy integrate the pointwise saddlepoint expansion of the density of 

(5.11) f{z) = -^exp{-h{ez)}!^l+^^p4{z) - ^fisizfj +en{z)y 
In (5.11) en{x) is a 0{n~'^) (for fixed z) function of z, and 

Mz)=r-'K^z\G^y^ PA{z)=T-'Kf{ez) 

are the standardized skewness and kurtosis of Z under the tilted (conjugate) 
density. The square-bracketed term in (5.11) is of order 0{n~^) for each fixed 
z. Hence, 

/)o,n(n)=P{Z(t)>n} 

(5.12) 

fn{z)dx+ / fn{z)[\pA{z) - ^^^{zf)dz, 



where fn{z) = eici>{—Iz{9z)} /W'^t^Tz] is assumed integrable. So provided 
(5.13) swp\\pa{z) - ^P3{zf + en{z)\ = o(n"°) 

z>u 

as n ^ oo, the first term in (5.12) will give an approximation with smaller 
relative error than pk,niu) or 'p^^{u) (assuming numerical integration does 
not compound the error). The tradeoff is an increase in computation, al- 
though solving the saddlepoint equation everywhere can be avoided by a 
judicious change of variables ([15], Section 6.2). It is well known (e.g., [15], 
Remark 3.2) that gains can be achieved when fn{z) is normalized to be- 

come a density. Let /„ {z) = /n(2)/||/n||i; a good approximation to po,n{u) 
is given by 

A A f°° 

P0,n (u) = / /„ (z) dz. 



z=u 
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For p-value calculation the order 0(n~") approximations />o,n and 
will often be sufficiently accurate. If the field is not too smooth and/or the 
search region is large in volume, higher terms begin to dominate (1.3) and the 
Gaussian density ^ may in fact be adequate. For a more complete account 
of tail probability approximations, see the references already mentioned in 
this subsection. 

Corollary 5 and Proposition 6 suggest that, asymptotically, tilting the EC 
densities is worthwhile provided u grows with n. In statistical applications, 
"li — > oo" corresponds to "large u," which is a crucial if unstated assump- 
tion when thresholding random field data via the expected EC heuristic. 
However, the tilted densities are not guaranteed to be any more accurate 
than the Gaussian ones in the small threshold scenario q = 1/2, at least as 
n — > oo. The performance of the approximation in finite samples is another 
matter altogether. These issues and others will be addressed in Sections 6 
and 7. 



6. Example: A scaled random field. It is instructive to consider a ran- 
dom field for which both the tilted and exact EC densities are available ex- 
plicitly; the Xn random field is one such case. Let Xi, X2, . . . , be i.i.d. cen- 
tred, unit variance Gaussian random fields on T and put Y{t) = Yl'i=i ^i(^)^- 
It has been established (by two distinct techniques [31, 38]) that the EC den- 
sities of Y are given by 



(A: - l)!n("-'=)/2e-"/2 
(2^)'=/2r(n/2)2"/2-i 

L(fe-1)/2J fc-l-2j 

k-l-2j-l 



(6.1) X E E ' 



j=o 1=0 



-^{n>k-2j-l} 



/!j!2i 
for A: > 1. 

We normalize by forming Z(t) = n~^^'^(Y{t) — n), so that Z is centered, 
approximately Gaussian, and = Var[Z(i)] = 2. Note that the G.eld is 
"twice as rough" as each Xi ; if the component Gaussian fields are isotropic 
with second spectral moment A^;, then (see [1], Chapter 7) Var(Z) = 
n~^Var(y) = 4Aa;lAr. To insure that A = 1 for the normalized process, we 

therefore assume that the Xj's are isotropic with Ax = Var(Xj) = In 1^- 
With that covariance structure and for a square parameter space T, we 
have 

E(^({t : Z{t) > u}) = Etf{{t : Y{t) > ^fRu -h n}) 
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k 
k 

Thus, the proper comparison is between pk^n 

(u) and 2 ^/O^ + for 

k = 0, . . . , N [this is a minor notational discrepancy; (6.1) is in fact 2^p^^{u) 
by the conventions of this article]. We have 

Kz{9) = -^e-'^\og( 1-2-^- 



K'z{0) 



29 



1-26/^' 



KUe) = 2{l-2-"^^ ' 



It follows that 



Vnu ^ , ' , ^/nu n , u , 



r^=2(l + 4=y; rX 



n/ y/2 

Note that the Hermite portion of the tilted density is identical to that of the 
Gaussian density in this case. Plugging these into (4.5), we have, for > 1, 

(6-2) PkA^) = (,^^^(.+i);4/2 H,_,{u/V2). 

Collecting all of the constants in (6.1), again for /c > 1, 

Pk,n{u)=2-''pl^{y/^U + n) 

(6.3) 

_ e-"v^/^(l + K/V?I)"/^-^ p 

(27r)(^+i)/22^/2 ^k-iA'^), 

where (assuming n > A^) Rk^n is the following degree k polynomial: 



^•"^ ' r n/2 \2) 



r(n/2) 

L^J^^Y n-l \ (-l)^+^+V+;-^ 
^ \k-2j-l) l\j\2j 



i=o /=o 

X ( 1 + 



X ^j+«+(fc+l)/2 
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Also, note that 

^J^^r^r^{n^i)l2 r(ra/2 + l) _ n 
r(n/2) yi) r(n/2) ~ 2 

as n — > oo. The corresponding Gaussian densities are 

(6-4) pIM = (^^^( fc>l. 

For the A; = EC density, we have 

Po,n{u) = Po,n(.^^ + n) = F{Y > y/nu + n}= gn{y) dy, 

' J y/nu+n 

where Y = Y{t) ~ xl and c/„,(y) = r(n/2)-i2-'^/2y'^/2-ie-2'/2 the density 
of Y. It is a weU-known fact [13] that the first-order saddlepoint approxima- 
tion Qn to gn is cxact up to a constant, so that 

(6.5) 9niy)=9n{y)- 

A 

To obtain a similar result for po,n('u) and Po.n (^^); we shall use the following 
lemma, which is easily proved. 

Lemma 7. Let Z = (Tn^(Y — be a random variable with density 
fn{z) = angni'^n^ + Pn) , where gn{y) is the density of Y . Then Mz E M, 

A A 

/„ (z) = Gn 9n {cTnZ + Pn)- 



Let fniz) denote the density of Z{t). Putting (T„ = y/n and fj,n = n in the 
lemma and using (6.5), we have 



PQ,n{Vnu + n)= 9n (y) dy 



n 



fAn~'/'{y-n))dy 

J y/nu+n 

oo A 

fn (z) dz 

= P0,n (u)- 

So the integrated, scaled tilted zeroth density is exact in this case, and this is 
a rare instance where normalization of the saddlepoint approximation does 
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not increase computation. The ordinary tilted and Gaussian A; = densities 
are given, respectively, by 



As a term-by-term comparison of (6.2) and (6.3) is onerous, we compare 
Hf^{u/^/2) and Rk,n{u) for k = 0,1,2, the cases of interest in a 2D or 3D 
problem. We have 



Ho{u/V2) = 1, i?o,n(u) = c„ Jl + ^ = 1 + 0(n-°); 



n 



„ , s U ( U 1 1\ U , ^, r.s-. 

n2 



2 lj(l + 0(n-")); 

where c„ = ^/2^e-"/2r(n/2)-i(n/2)("-i)/2 = 1 + 0{n^^). Hence, the tilted 
EC densities for A; = 1,2,3 have a relative error 0(n~"). Note that the 
Gaussianity of the Xj's in this example leads to K[Z'^Zj] = 0, which im- 
plies that, for all A; > 1, a relative error no worse than 0{n^^^^^'^'^^'^~^'^^) 
must be achieved by pk^n — see (4.7). It is easily verified that luidu) = u'^/4: — 
v?/{Q^/n) + 0(n^~^"), so that en[/Ofc^„('"), Pfc,n(w)] = 1 - 3a when a > 1/3. 

The behavior of these approximations for A; < 3 is explored further in 
Figure 2. In each column we consider a different asymptotic regime: the 
most general one in which a = 1/3 and both n and u grow; the classical one 
in which a = 1/2 so that u is fixed while n grows; and the less common one in 
which n is fixed while u grows (corresponding to a = — oo). The latter regime 
is meaningful for statistical applications since the sample size (ra) is large 
but finite and the relative error in '&'^{Au) as a p- value approximation is 
conjectured to fall off exponentially as n — > oo (see Section 4.4). The theory 
of Sections 4.1 and 5 does not directly apply to this scenario. A fourth 
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Fig. 2. iJC densities for the Xn random field: exact (pk,n, " — "); tilted (pk,n, "- -"); and 
Gaussian (pi ^, "■■■")■ Ratios r of approximate to exact densities are plotted for three 
regimes: u ~ cn"''"' oo (a — 1/3, 1st col.); u fixed and n ^ oo (a — 1/2, 2nd col.); n 

A ^ _ 

fixed and it — > cxa (a = —oo, 3rd col.). For k = 0, po.n =Po,n, and "-• -" denotes „. 
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regime, where /c — > oo, might also be of interest; empirical evidence for the 
Xn field suggests that the tilted densities perform poorly under that scheme. 

For k = 0,l in this example, the tilted densities offer a clear improvement 
over the Gaussian approximation even for small n and/or u. For k = 2,3, 
a larger sample size or threshold is needed before the difference emerges, 
but it is still unequivocal when it — > oo. When a = 1/2 for k = 2 or 3, 
the tilted and Gaussian approximations perform similarly, as predicted by 
Corollary 5. In the third regime where u — > oo while n remains fixed, the 
Gaussian densities are quickly diverging from the correct densities; however, 
they are more accurate than the tilted densities for k = 2,3 in the range 
2 < It < 4. Note that the chosen sample size n = 500 was quite large in that 
regime. When k = the naive tilted density po.n(^) appears to outperform 
the Lugannani-Rice estimate's leading term 7jo„('u), but this difference is 
small and may be specific to the Xn field. The Gaussian assumption gener- 
ally leads to under estimation of the expected EC in this example. In other 
words, using a Gaussian assumption for statistical thresholding of Xn fields 
is anti-conservative. 

7. Application: Subtracted "bubbles" images. A moment's thought re- 
veals that the form (1.1) taken for the Z process is far more restrictive than 
is necessary. Indeed, the crucial property of Z = (Z, Z' , Z') appearing in the 
proof of Theorem 4 is the rate of decay of its (joint) cumulants; namely 

(7.1) E[Z(i)]=0, Var[Z(t)] =cj2; 
and 

(7.2) A:,([Z(i)r)=0(n-(^"2)/2)^ 7 G n(^+^)(^+2)/M7| = z., > 2. 

One might therefore expect Theorem 4 to apply to a variety of stationary 
random fields satisfying (7.1) and (7.2) in addition to the regularity con- 
ditions in the Appendix. The following is a slightly more general construc- 
tion than (1.1): let {Wi(t)}i be independent but perhaps nonidentically- 
distributed random fields with finite cumulants of all orders. Assume that 
the random field of interest takes the form 

1 " 

(7.3) Z(t) = ^5]a,W^,(t) 

for finite constants (0^)^^^. This framework covers the case of a two-sample 
problem (e.g., a difference in means between two groups), as well as the 
following application. 

The bubbles methodology of Gosselin and Schyns has been described in 
detail elsewhere [16]. Typically, the experimenters' goal is to determine which 
parts of an image are most important in a visual discrimination task. The 



EULER CHARACTERISTICS OF CENTRAL LIMIT FIELDS 27 



(a] Happy (b) Fsgrful 




50 100 S50 200 250 50 100 150 200 260 



(€) Biibb^es (dj Stimulus 




50 100 150 £00 250 50 100 150 200 250 



Fig. 3. Bubbles experiment. The subject is asked to discriminate between the happy (a) 
and fearful (b) faces on presentation of the stimulus (d) which is one of the two faces, 
here the fearful face, partially revealed by random "bubbles" (c). The 227^ search region T 
is inside the black frame in (a) and (b) . 

data analyzed later in this section come from [3] . On each trial of the exper- 
iment, a subject was shown a 256^ image of a face that was either fearful or 
happy. The images were masked apart from a random number of localized 
regions or "bubbles" that revealed only selected parts of the face; see Fig- 
ure 3. The subject was then asked whether the partially revealed face was 
fearful or happy, and the trial was repeated about 3,000 times on each of 10 
subjects. 

Each masked image was generated as follows. First a range of scales was 
chosen, asj = 3 x 2-' , j = 1, . . . , 5. Then the original image Xq [Figure 3(a) 
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or (b)] was smoothed by an isotropic Gaussian filter with standard devia- 
tion a"sj/2 to produce images Tj. The smoothed images were differenced to 
produce images Dj = Ij-i — Ij that reveal image features at five different 
scales. Differenced image Dj was then multiplied by a mask consisting of 
the sum of a random number of isotropic Gaussian "bubbles," each with 
standard deviation crsj- The bubble center were chosen at random from the 
256^ pixels, that is, according to a Poisson process. The number of bubbles 
for each scale was a multinomial random variable with probabilities inversely 
proportional to bubble area (tx cr^"^), so that equal areas were revealed by 
each of the five bubble sizes, on average. The sum of all the bubbles over all 
five scales is shown in Figure 3(c) for one of the trials. The sum of all the 
bubbles times the differenced images is shown in Figure 3(d). On the basis 
of Figure 3(d) the subject must decide if the face is happy or fearful (fearful 
is the correct answer in this case). The total number of bubbles was chosen 
dynamically to maintain a roughly 0.75 probability of correctly identifying 
the face. 

7.1. Test statistic, cumulants and tilted densities. For the present anal- 
ysis we consider a single-subject run. We assume that the grid is square, 
and that the bubble standard deviation as (or equivalently, the bubble full 
width at half max F = \/81og2cj_B) is fixed, as is the number m of bubbles 
per image. We denote by pc and pi = 1 — pc the proportions of correctly 
and incorrectly classified images. Let Yip denote the number of bubbles in 
image i which are centred at pixel p, with i = 1, . . . , n and p = l,. . . ,P. Then 
Yip r^V{m/P) is a Poisson random variable, and the Yip are approximately 
independent (exactly when they come from different images). The test statis- 
tic field for detecting significant regions is formed by signing and summing 
the bubbles, with sign depending on whether the image was correctly or in- 
correctly classified. Thus, the unnormalized test statistic at pixel t is given by 



n p 



(7.4) 



Z{t) = n- 



i=l p=l 



where 




image i correct, 
image i incorrect 



and 




B 



The statistic is normalized by forming 



Z(t) = Z{t)l^\lzxZ{t). 
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Inference can be either conditional on the observed vector a or uncondi- 
tional, in which case the are i.i.d. scaled Rademacher random variables 
and PC is taken as the (assumed known) population proportion. One may 
apply the theory of the present article under either analysis; we have chosen 
the conditional approach since it simplifies cumulant calculations. Thus, pc 
and pj are observed proportions. 

The most natural global null hypothesis is 

Hq : Oj is independent of (lip)^=i Vi < n, 

which implies 

HQ:E{Z{t)\a} = ^teT, 

since EH,{Z{t)\a}cxJ:pbpit)J:ia^lKHo{Y^p\a} = ^(Ep bp(t))(Ei a*) = 0. Ig- 
noring mild boundary effects, note that the null distribution of Z(t) is 
unchanged if we translate or rotate the coordinate system; thus, Z{t) is 
(strictly) isotropic. By (7.4), Z is of the form n'^^'^J^i^^i^i where the Wj's 
are i.i.d. random fields. The distribution of Wi{t) = Ep^p(0^p is very non- 
Gaussian, with a high probability of being nearly zero (if there are no bubble 
center near t), and a small probability of being large. The distributions be- 
come more Gaussian as m or F increase. 

An important practical note is that most users of this methodology choose 
to "clamp" the summed bubbles at a maximum height equal to the peak 
height of a single bubble. In particular, two bubbles centered at the same 
pixel, when added, have the same peak height as a single bubble. This 
is how the mask is actually created, and it is felt that this is how the 
stimulus is perceived. Clamping effectively replaces (7.4) with 
n~^/'^J27=i^i^^^{^p=ibpi^)Yip^^}- Doing so not only complicates cumu- 
lant derivations, but, more importantly, degrades the "smoothness" of the 
random field. When the number of bubbles is small compared to the number 
of pixels and their width is not too large (as in the data analyzed below), the 
loss of smoothness due to clamping is negligible. Our choice of an undamped 
analysis is always valid, but perhaps less powerful at detecting facial features 
if indeed the probability of response is related to the clamped rather than 
the undamped bubble mask. 

Without any extra effort we generalize the cumulant computations to N 
dimensions, although currently in bubbles experiments N = 2. It can be 
shown that the jth null conditional cumulant of Z{t) is approximately 

N/2 pi-i ^ (-i)ip}-J ^ / 41og2 \^/^ 1 

which only depends on pc and the total number B = nm/{PF~^) of bub- 
bles per A^-dimensional "resel" (volume element). These have been used to 
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compute the tilted variance of Z, 

In practice (for typical values of pc and i?), j < 20 seems to be enough; 
beyond this the cumulants are extremely close to zero. The tilting parameter 
6u is found by solving K'^{9) = u with a Newton-Raphson algorithm started 
from the Gaussian solution 9 = u; five iterations appear to suffice in most 
cases. An approach like this one was advocated by RS. 

One can also show that the second spectral moment of Z under Hq is 
given (approximately) by 

(7.6) Var{Z{t)\a}^-^lN = XlN. 
The p- value approximation when = 2 is thus 

(7.7) pjsup Z{t) > n| « E^{A^) = po{u) + 2{PX)^/^pi{u) + P^Mu), 

where the EC densities are 

Mu) = - Iu{0u)}'f{rJu), 

(7.8) ^i(n) = exp{-/„(^„)}(27rT„)-S 
P2{u) = exp{-/„(^„)}(27r)-3/V-ie"„. 

By contrast, the Gaussian theory EC densities are 
p2{u)=^{u), 

(7.9) pUu)=exp{-u'/2}i27r)-\ 
p](n)=exp{-t.V2}(2vr)-3/2^, 

with p- value approximation 

(7.10) p2iu)+2iP\)y'pJ{u)+PXpl{u). 

When P is large, the third terms in (7.7) and (7.10) dominate. Typical values 
for the parameters are P = 256^, m = 16.5, F = 14.1, ?i = 3000, pc = 0.75; 
setting u = 3.965 [the p = 0.05 threshold using (7.7)] then gives 

Po{u) = 1.4 X 10-^ 2{PX)'^/^pi{u) = 0.00171, PXhH = 0.0484. 

In view of this, the ratio of the tilted to Gaussian p-value approximations, 

(7.11) r = exp\ — - Iu{Ou) > , 

P2W 12 } TuU 
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only depends on u, B and pc- We also consider the approximation suggested 
by extension (4.18) of RS, namely, 

(7.12) Mu) + 2(PA„)i/2pi(n) + PKMu), 

where A„ denotes the tilted second spectral moment of Z. It can be demon- 
strated that Au ~ r^A. To compare this approximation to the Gaussian the- 
ory, one can again form the p- value ratio 

(7.13) rRs ^ ^ w r^r. 

These ratios are plotted for u = 3.5 and a range of B and pc in Figure 4. 
Note that in this application the approximation (7.7) deviates less extremely 
from the Gaussian one than does the heuristic of RS. This can be understood 
by noting that, when a = 1/2, 

(7.14) ^2^1 + ^3^^ + 0(^-1), 

where 6*^ ~ > and sgnjKs} = sgn{| — pc}- At the same time, r > 1 44> 
Pc < 1/2. The 0{n~^) term in (7.14) can be positive or negative but is of 
smaller order than the skewness term. It is also worth noting that unlike in 
the Xn example and the fields considered by RS, the Gaussian approximation 
can be conservative (i.e., lead to larger values) in bubbles experiments 
when typical parameter values {pc > 1/2) are used. 



7.2. Simulation. A simulation was devised to compare (7.7), (7.10) and 
(7.12) with the true tail probability (1.2) of sup^ Z[t) above a threshold un- 
der Hq. In each of 10^ Monte Carlo iterations, we generated 3000 bubbles 
images on a P = 208 x 208 square with corresponding classifications, using 
Pc = 0.75 and bubble width F = 24. For each iteration we then calculated a 
sequence of test statistic fields Zn{t) for n = 100, 200, . . . , 3000. Since kernel 
smoothing with Fourier methods was used to speed up computation, the data 
were generated first on a larger 256 x 256 square before discarding the outer 
edges to eliminate periodicity. The number of bubbles per image was fixed so 
that an average of m = 20 would fall within the search region; this has a neg- 
ligible effect on the cumulants (7.5). The mixed skewness terms E[Z(f)^Zj(f)] 
can be shown to be (approximately) proportional to Yl,p{ti ~ Pi)^p{tY ^ and 
hence, are very close to zero for pixels far from the boundary (recall that 
they were exactly zero in the example). We therefore expect the tilted 
approximations to be within the order 0(n~°) of the true p-value. Figure 5 
displays the simulation results under the three asymptotic regimes. In ad- 
dition to using both A and A„, we also tried u/a in the place of TuOu as 
the argument in the Hermite polynomials. This resulted in third and fourth 
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Fig. 4. Ratios of tilted to Gaussian-approximated expected Euler Characteristic of the 
excursion of a test statistic from a typical bubbles experiment above the threshold u = 3.5, 
using both the untilted (\) and the tilted (Xu) second spectral moment. The ratios are 
plotted for a reasonable range of B (the total number of bubbles per resel) and for several 
values of PC (the proportion of correctly classified images). When pc = 1/2 (red curve), 
the skewness of Z is zero, and the saddlepoint correction is less pronounced. Assuming 
that the tilted expressions are closer to the true p-value, r > 1 signifies that the Gaus- 
sian approximation is anti-conservative, while r < 1 means that the Gaussian procedure is 
conservative. Note that pc ~ 0.75 is commonly targeted by the originators of the bubbles 
paradigm, making the Gaussian theory valid though less powerful than the tilted inference. 

variations on the tilted p-value approximation. In all instances the tilted ex- 
pected EC's were closer to the observed p-value than the Gaussian expected 
EC. In particular, the expansion which uses the tilted spectral moment 
and the Gaussian Hermite argument u/a (" — ," orange) performed best, al- 
though differences between tilted approximations were small. The observed 
directional trend between the four tilted expressions is likely specific to this 
application. We note that while feasible for parameter spaces of this size, 
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deriving p-values by simulation for bubbles-type data is time-consuming, 
and quickly becomes impracticable as the dimension N increases beyond 2. 
The benefits of having a closed-form p- value approximation are evident. 

7.3. Data analysis. We analyzed a portion of the data from [3] described 
at the beginning of this section (control group). One subject (#9) was pre- 
sented with n = 3072 images of size 256 x 256 pixels, which he/she classified 
as either happy or fearful, at a rate of pc = 0.76. The search region was 
trimmed to the inner P = 227^ pixels to avoid Fourier edge effects. In agree- 
ment with the theory in this article, we have limited our analysis to a single 
bubble width, F = 18.8 pixels or as = 8 (the subject was in fact shown bub- 
bles of width F = 28.3, i.e., j = 2, asj = 12; but we chose a narrower filter 
when smoothing the test statistic in order to totally eliminate the "clamp- 
ing" issue). The mean number of bubbles per image at this scale and within 
the search region was m = 3.0. Our chosen single bubble width was the 
second-narrowest class (j = 2); an example of the three bubbles is clearly 
visible in Figure 3(c). These parameter values correspond to B = 63.6 bub- 
bles/resel. Thresholded images are displayed in Figure 6. There is strong 
signal in the right eye region for this subject. The global maximum of Z{t) 
exceeded 5.8, and hence, all thresholding methods identified that peak. At 
the p = 0.05 level the Gaussian threshold picked up only the right eye, while 
the tilted threshold found additional activation in the upper lip. At p = 0.1 
significance the tilted expected EC identified a third signal in the left eye 
which went undetected by the Gaussian theory. Note that the Gaussian ap- 
proximation is conservative in this case since Z(t) exhibits negative skewness 
(PC > 1/2). 

8. Discussion and future work. We have provided a rigorous theoretical 
justification for the use of saddlepoint methods in thresholding asymptotically- 
Gaussian random fields, first proposed by Rabinowitz and Siegmund [25]. In 
particular, we have extended their "expected number of maxima" heuristic 
to an expected Euler Characteristic approximation for locally isotropic Cen- 
tral Limit- type fields, and derived relative errors for the EC densities when 
the threshold is allowed to grow. We have also compared these formulae to 
a Gaussian approximation, both analytically and in numerical and real data 
examples involving a Xn field and a linear combination of Poisson fields. 

As did RS (who only considered a Poisson process) vis-a-vis expected 
local maxima, we found substantial differences between the Gaussian and 
the various tilted approximations (which were uniformly more accurate in 
a particular case — see Figure 2) to the expected EC. It was shown in Sec- 
tion 5 that in a large-threshold setting the saddlepoint correction reduces 
asymptotic relative error. In contrast to the examples chosen by RS, we have 
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Fig. 5. Ratios of expected EC approximations to true (empirical) p-values for simulated 
bubbles data, under the three asymptotic regimes of Section 7. Bubbles parameter values 
were P = 208^, F = 24, m = 20, pc = 0.75, and n ranging from 100 to 3000. The ratio scale 
was chosen so that absolute deviations from the line r = 1 correspond to relative error. The 
Gaussian expected EC is given by (7.10). The tilted approximations use either the untilted 
(X) or tilted (\u) second spectral moment, and either the tilted (tuOu) or untilted (u/a) 
Hermite argument. Monte Carlo error bars (±1 standard deviation) are for the Gaussian 
curve ("■ ■ ■ "), but are attached to the line r — 1 to aid interpretation. They are computed 
pointwise, using \/sr(jp/pemp) ~ 10^^p^(l — Pemp)/Pemp) where pemp is an empirical p-value 
and p one of the EC approximations. Error bars for the other curves are slightly narrower. 
The threshold u — 3.5 was chosen in (b) in order to have empirical probabilities near 0.05. 
The sample size n = 300 was chosen in (c) m order to have the total bubbles/resel close to 
the value B = 63.6 observed in the data. In all cases the tilted approximations outperform 
the Gaussian one, including in (c) where n is fixed, u grows and all curves appear to 
diverge from unity. 
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Fig. 6. Thresholded bubbles data for subject #5, with the search region framed. Test 
statistic (a), with excursions above: (b) u = 2 (Gaussian uncorrected p ~ 0.05^; (c) u — 3.92 
(Gaussian pr^ 0.05); (d) u = 3.76 (tilted p ^ 0.05); (e) u = 3.72 (Gaussian p r=i 0.1); (f) 
M = 3.58 (tilted p ^ 0.1). P-values approximated by KEC under the Gaussian (7.10) and 
tilted theory with argument t^Ou and spectral moment A„; see (7.12). Inference using (7.7) 
is similar and omitted. There is strong signal in the right eye region. At p = 0.05, the 
Gaussian theory identifies only the right eye (c) ,• tilting finds activation on the upper lip 
(d). Atp = 0.1, tilting detects signal above the left eye (f); the Gaussian approximation does 
not (e). Tilting reduces thresholds because pc > 1/2, so Z(t) exhibits negative skewness. 
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demonstrated via simulation that the Gaussian approximation can be con- 
servative in some instances rather than always underestimating the p- value; 
see Section 8 and Figures 4-5. We have shown that the important correction 
for non-Gaussianity occurs in the exponential portion of the EC densities. 
Additional corrections to the variance of the first derivative field and to the 
argument of the Hermite polynomial are of smaller order and may not be 
needed unless the expected EC over the full range of u is required. The 
examples considered indicate that tilting the spectral moment has a greater 
impact than tilting the Hermite argument. Our results support RS's recom- 
mendation to employ a tilted expansion when the additional computation is 
feasible; we add that the benefit is greatest if the field has nonzero skewness 
and/or if u — > oo (is large) along with n. 

In addition to the bubbles test statistic field presented above, we have 
begun to apply the theory of this article to lesion density maps from neu- 
roimaging [10], which can be non-Gaussian in a similar fashion. Handling 
this and further applications will require generalizing these results to certain 
nonisotropic and nonstationary CLRF's; this work is already underway. As 
discussed in Section 4.4, it would also be desirable to rigorously extend the 
proof in [33] of the EC heuristic to the asymptotically Gaussian case. In a 
forthcoming paper we shall explore a geometric interpretation to the tilted 
EC densities in the spirit of [31]. This will tie into similar investigations of 
"CL-related" processes, including asymptotically t and F random fields. 

APPENDIX: REGULARITY CONDITIONS FOR THEOREM 4 

We assume that the random vectors W* = (VF*, vec[<TW^*|fc_i -|- 
W^'I/a]')' are nonlattice. (While saddlepoint expansions have been estab- 
lished for lattice distributions, the resulting random fields could scarcely 
be imagined to be regular in the manner described below.) There are sim- 
ple conditions which validate the tilted expansion (2.6) of the joint density 
f{z,z,m) of {Z{t), Z\k(ty , M{ty) [which exists, and is bounded and contin- 
uous, under said conditions]. For example, it is sufficient that: (i) W* come 
from an exponential family indexed by € O an open, convex subset of M'^*^ 
[which appears tacitly in (2.2)]; and (ii) W have characteristic function Tpg 
for which |'i/;5i(s)|'' is integrable for some v = 1^(9) > 1. For the latter, com- 
pare [4], Theorem 6.5, and their condition (c'). An extension of Theorem 
4 to the case where the Wi's are nonidentically distributed could also be 
conceived; for regularity conditions the reader is referred to Theorem 6.6 of 
the same monograph. 

The exponential family assumption guarantees that the W* have finite 
joint cumulants of all orders, but to bound the error we shall assume, in 
particular, that for every v <A the cgf of W* belongs to C^iV) for some 
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neighborhood V containing 0. We shall also require that for some finite 
constants H > 0, {bi)i>2 and all s>2, 

(A.l) E\\W\\' <'^b'^H'-^. 

As for the validity of (1-5), conditions are laid out in [1] and AT as follows: 

A T G is compact and dT has zero Lebesgue measure. 

B Z is almost surely suitably regular on T at the level u, namely: 

1. Z has continuous partial derivatives up to second order in an open 
neighborhood of T; 

2. there is no t G T such that Z{t) — u = Zi{t) = ■ ■ ■ = Z}^{t) = 0; 

3. there is no t G dT and permutation p of {1, . . . , /c} such that 

Z{t)-u = Zp,{t) = --- = Zp^_^{t)=Q; 

4. there is no t G T and permutation p of {1, . . . , /c} such that 

Z{t)-u = Zp^{t) = --- = Zp^_^ (t) = det(Zp,p^. {t))ij<k-i = 0. 

C The second-order derivatives of Z have finite variances. 

D The joint density of (Z, is continuous in each argument. 

E The conditional density of {Z, Z\'^_^) given (Z^, is bounded above. 

F The moduli of continuity LOi and tOij of the various derivatives satisfy 

pimax[a;i(?7),a;ij(77)] > e\ = oij]^) 

as r/ — > for every e > 0. 

Note that under the oft-made Gaussian assumption, many of these condi- 
tions are trivially satisfied. 
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